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MULTI-ADAPTIVE GALERKIN METHODS FOR ODES I* 

ANDERS LOGGt 

Abstract. We present multi- adaptive versions of the standard continuous and discontinuous 
Galerkin methods for ODEs. Taking adaptivity one step further, we allow for individual time- 
steps, order and quadrature, so that in particular each individual component has its own time-step 
sequence. This paper contains a description of the methods, an analysis of their basic properties, and 
a posteriori error analysis. In the accompanying paper [A. Logg, SIAM J. Sci. Comput., submitted], 
we present adaptive algorithms for time-stepping and global error control based on the results of the 
current paper. 
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1. Introduction. In this paper, we present multi-adaptive Galerkin methods 
for initial value problems for systems of ODEs of the form 

f u(t) = f(u(t),t), ie(0,T], 
[ ' 1 u(0) - u , 

where u : [0, T] -» R N , f : R N x (0,T] — > R N is a given bounded function that is 
Lipschitz-continuous in u, uo £ M. N is a given initial condition, and T > is a given 
hnal time. We use the term multi-adaptivity to describe methods with individual 
time-stepping for the different components Ui(t) of the solution vector u(t) = (ui(t)), 
including (i) time-step length, (ii) order, and (iii) quadrature, all chosen adaptively 
in a computational feedback process. In the companion paper [29], we apply the 
multi-adaptive methods to a variety of problems to illustrate the potential of multi- 
adaptivity. 

The ODE (|1.1[) models a very large class of problems, covering many areas of 
applications. Often different solution components have different time-scales and thus 
ask for individual time-steps. A prime example to be studied in detail below is our 
own solar system, where the moon orbits around Earth once every month, whereas 
the period of Pluto is 250 years. In numerical simulations of the solar system, the 
time-steps needed to track the orbit of the moon accurately are thus much less than 
those required for Pluto, the difference in time-scales being roughly a factor 3,000. 

Surprisingly, individual time-stepping for ODEs has received little attention in the 
large literature on numerical methods for ODEs; see, e.g., [H [21] [22j [3] [34] . For specific 
applications, such as the n-body problem, methods with individual time-stepping have 
been used — see, e.g., [3TJ H [5] or [25] — but a general methodology has been lacking. 
Our aim is to fill this gap. For time-dependent PDEs, in particular for conservation 
laws of the type ii + f(u) x = 0, attempts have been made to construct methods with 
individual (locally varying in space) time-steps. Flaherty et al. [20] have constructed 
a method based on the discontinuous Galerkin method combined with local forward 
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Euler time-stepping. A similar approach is taken in [BJ, where a method based on 
the original work by Osher and Sanders [33] is presented for conservation laws in one 
and two space dimensions. Typically the time-steps used are based on local CFL 
conditions rather than error estimates for the global error and the methods are low 
order in time (meaning < 2). We believe that our work on multi-adaptive Galerkin 
methods (including error estimation and arbitrary order methods) presents a general 
methodology to individual time-stepping, which will result in efficient integrators also 
for time-dependent PDEs. 

The methods presented in this paper fall within the general framework of adap- 
tive Galerkin methods based on piecewise polynomial approximation (finite element 
methods) for differential equations, including the continuous Galerkin method cG(q) 
of order 2q, and the discontinuous Galerkin method dG(g) of order 2q + 1; more 
precisely, we extend the cG(q) and dG(g) methods to their multi-adaptive analogues 
mcG(g) and mdG(q). Earlier work on adaptive error control for the cG(g) and dG(q) 
methods include [16j [24J [TH [TT1 [19] . The techniques for error analysis used in these 
references, developed by Johnson and coworkers (see, e.g., [TTJ [12l [101 H3J HH [15] , and 
[8] in particular) naturally carries over to the multi-adaptive methods. 

The outline of the paper is as follows: In section[2]we summarize the key features 
of the multi-adaptive methods, and in section [3J we discuss the benefits of the new 
methods in comparison to standard ODE codes. We then motivate and present the 
formulation of the multi-adaptive methods mcG(g) and mdG(q) in section [4] Basic 
properties of these methods, such as order, energy conservation, and monotonicity, are 
discussed in section[5J In the major part of this paper, section[BJ we derive a posteriori 
error estimates for the two methods based on duality arguments, including Galerkin 
errors, numerical errors, and quadrature errors. We also prove an a posteriori error 
estimate for stability factors computed from approximate dual solutions. 

2. Key features. We summarize the key features of our work on the mcG(g) 
and mdG(q) methods as follows. 

2.1. Individual time-steps and order. To discretize (ll.ip . we introduce for 
each component, i = 1,,..,N, a partition of the time-interval (0,T] into Mi subin- 
tervals, iy = (ti j—i, ty], j = 1, . . . , Mi, and we seek an approximate solution U(t) 

I such that Ui(t) is a polynomial of degree qij on every local interval iy. Each 
individual component Ui(t) thus has its own sequence of time-steps, {fcy}*LV The 
entire collection of individual time- intervals {iy} may be organized into a sequence 
of time-slabs, collecting the time-intervals between certain synchronised time-levels 
common to all components, as illustrated in Figure |2~T1 

2.2. Global error control. Our goal is to compute an approximation U(T) of 
the exact solution u(T) at final time T within a given tolerance TOL > 0, using a 
minimal amount of computational work. This goal includes an aspect of reliability 
(the error should be less than the tolerance) and an aspect of efficiency (minimal 
computational work). To measure the error we choose a norm, such as the Euclidean 
norm || ■ || on R^, or more generally some other quantity of interest (see |32|). 

The mathematical basis of global error control in || ■ || for mcG(g) is an error 
representation of the form 



where R = (Ri) = R(U,t) = U(i) — f(U(t),t) is the residual vector of the approximate 



(2.1) 
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Fig. 2.1. Individual time- discretizations {or different components 



solution U(t), ip(t) is the solution of an associated linearized dual problem, and (•, •) 
is the M. N scalar product. 

Using the Galerkin orthogonality, the error representation can be converted into 
an error bound of the form 

N 

(2.2) \\U(T)-u(T)\\ <$>(T) ^Utr^iUM 

i=i 

where {Si(T)}^_ 1 are stability factors for the different components, depending on the 
dual solution <p(t), and where ki{t) = kij, <fc(i) — q^ for t £ Iij. The error bound may 
take different forms depending on how J Q (R, ip) dt is bounded in terms of R and <p. 

By solving the dual problem numerically, the individual stability factors Sj(T) 
may be determined approximately, and thus the right-hand side of (|2.2p may be 
evaluated. The adaptive algorithm seeks to satisfy the stopping criterion 

N 

(2.3) ^Si(T) mg ki(t)«W\Ri{U,t)\ < TOL, 
»=i — 

with maximal time-steps k = (ki(t)). 

2.3. Iterative methods. Both mcG(g) and mdG(g) give rise to systems of 
nonlinear algebraic equations, coupling the values of U (t) over each time-slab. Solving 
these systems with full Newton may be quite heavy, and we have instead successfully 
used diagonal Newton methods of more explicit nature. 

2.4. Implementation of higher-order methods. We have implemented 
mcG(g) and mdG(q) in CH — h for arbitrary q, which in practice means 2q < 50. The 
implementation, Tanganyika, is described in more detail in [29] and is publicly (GNU 
GPL) available for Linux/Unix [50] . 

2.5. Applications. We have applied mcG(g) and mdG(q) to a variety of prob- 
lems to illustrate their potential; see [29]. (See also [27] and [26].) In these applica- 
tions, including the Lorenz system, the solar system, and a number of time-dependent 
PDE problems, we demonstrate the use of individual time-steps, and for each system 
we solve the dual problem to collect extensive information about the problems stability 
features, which can be used for global error control. 
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3. Comparison with standard ODE codes. Standard ODE codes use time- 
steps which are variable in time but the same for all components, and the time-steps 
are adaptively chosen by keeping the "local error" below a given local error tolerance 
set by the user. The global error connects to the local error through an estimate, 
corresponding to (|2.2I) , of the form 

(3.1) {global error} < S max{local error}, 

where S is a stability factor. Standard codes do not compute S, which means that 
the connection between the global error and the local error is left to be determined 
by the clever user, typically by computing with a couple of different tolerances. 

Comparing the adaptive error control of standard ODE codes with the error con- 
trol presented in this paper and the accompanying paper |2S], an essential difference 
is thus the technique to estimate the global error: either by clever trial-and-error or, 
as we prefer, by solving the dual problem and computing the stability factors. Both 
approaches carry extra costs and what is best may be debated; see, e.g., [3 2) for a 
comparison. 

However, expanding the scope to multi-adaptivity with individual stability factors 
for the different components, trial-and-error becomes very difficult or impossible, and 
the methods for adaptive time-stepping and error control presented below based on 
solving the dual problem seem to bring clear advantages in efficiency and reliability. 

For a presentation of the traditional approach to error estimation in ODE codes, 
we refer to [2] , where the following rather pessimistic view is presented: Here we just 
note that a precise error bound is often unknown and not really needed. We take the 
opposite view: global error control is always needed and often possible to obtain at a 
reasonable cost. We hope that multi-adaptivity will bring new life to the discussion 
on efficient and reliable error control for ODEs. 

4. Multi-adaptive Galerkin. In this section we present the multi-adaptive 
Galerkin methods, mcG(g) and mdG(<7), based on the discretization presented in 
section 12.11 

4.1. The mcG(q) method. The mcG(g) method for (jl.ljl reads as follows: 
Find U eV with 17(0) = Uq, such that 

(4.1) f (U,v)dt= [ (f(U,-),v)dt WveW, 

Jo Jo 

where 

U9) V = {ve[C([0,T\)] N e7>»* (hj), 3 = !,-■•, M u i = 1,..., N}, 

1 ' ' W = {v. uil/y G V^-Hlii), 3 = 1. • • • ,M t , i = 1, . . . , N}, 

and where V q (I) denotes the linear space of polynomials of degree < q on I. The 
trial functions in V are thus continuous piecewise polynomials, locally of degree g^, 
and the test functions in W are discontinuous piecewise polynomials that are locally 
of degree qij — 1. 

Noting that the test functions are discontinuous, we can rewrite the global prob- 
lem (|4.1|) as a number of successive local problems for each component: For i = 
l,...,N,j = l,...,Mi, find Ui\ hj G V q ^ (hj) with given, such that 



(4.3) 



/ u iV dt= [ fi(u,-)vdt Vver q ^-\iij). 

J I J I 
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We notice the presence of the vector U(t) = (Ux(t), . . . , Ujy(t)) in the local prob- 
lem for Ui(t) on Iij. If thus component U ix (t) couples to component Ui 2 (t) through /, 
this means that in order to solve the local problem for component Ui 1 (t) we need to 
know the values of component Ui 2 (t) and vice versa. The solution is thus implicitly 
defined by (14.31) . Notice also that if we define the residual R of the approximate 
solution U as Ri(U,t) = Ui(t) — fi(U(t),t), we can rewrite (|4.3p as 



(4.4) 



/ Ri{U, -)v dt = VveV^-^Iij), 
J hi 



i.e., the residual is orthogonal to the test space on every local interval. We refer to 
this as the Galerkin orthogonality for the mcG(g) method. 

Making an ansatz for every component Ui(t) on every local interval Zy in terms 
of a nodal basis for V qii (Iij ) (see the appendix) , we can rewrite (|4.3[) as 

(4.5) £y m = £y + / w [ £ i] (th (t)) fi(U(t),t) dt, m = 1, 



J hi 

where {£ym}m=o are the n °dal degrees of freedom for Ui(t) on the interval Zy, 
{u>m }m=i <- V q ~ l (Q, 1) are corresponding polynomial weight functions, and ry maps 
Iij to (0,1]: Tij(t) = [t — t%j—x)/{t%j — k,j-i)- Here we assume that the solution 
is expressed in terms of a nodal basis with the end-points included, so that by the 
continuity requirement £yo = Ci,j-i,qt <_i< 

Finally, evaluating the integral in (|4.5p using nodal quadrature, we obtain a fully 
discrete scheme in the form of an implicit Runge-Kutta method: For i = l,...,N, 
j = 1, . . • , Mi, find {£,ijm}m J =cn with £yo given by the continuity requirement, such 
that 

(4.6) £y m = + % £ "fe 1 /«(C/(^ 1 (4 9 « ] )) I Ty 1 (4«])), m = 1,.. . ,qy, 

n=0 

for certain weights {iOm n } and certain nodal points {s^} (see the appendix). 

4.2. The mdG(q) method. The mdG(g) method in local form, corresponding 
to (|3~3"j) . reads as follows: For i = 1, . . . , N, j = 1, . . . , M u find Ui\ hj € P 9 « (Zy), such 
that 

(4.7) [U i ] i ^ 1 v(tf j _ 1 )+ f UiVdt= f fi(U,-)vdt WeP 9 «(Zy), 

where [•] denotes the jump, i.e., [u]y = u(iy) — w (*y) ; an d the initial condition is 
specified for i = 1, . . . , N, by ?7j(0 _ ) = u,(Q). On a global level, the trial and test 
spaces are given by 

(4.8) V = W = {v:v i \ Iij ■ V" (/,,:../ ! \/,. i = l,...,N}. 

In the same way as for the continuous method, we define the residual R of the ap- 
proximate solution U as Ri(U,t) = Ui(t) — fi(U(t),t), defined on the inner of every 
local interval Zy , and we rewrite (|4. 7[) in the form 



(4.9) [U i ] i , j -iv(t+) + / Ri{U, > dt = Vve (Zy ). 
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We refer to this as the Galerkin orthogonality for the mdG(g) method. Notice that 
this is similar to (|4.4j) if we extend the integral in (|4.4I) to include the left end-point of 
the interval Jy . (The derivative of the discontinuous solution is a Dirac delta function 
at the end-point.) 

Making an ansatz for the solution in terms of some nodal basis, we get, as for the 
continuous method, the following explicit version of (|4.7j) on every local interval: 



(4.10) 6im = £y + / tC* J (ry (*)) fi(U(t), t)dt, m = 0, . . . , q ih 



(4.11) Z iSm = + h 3 £ MU(Tr/(s^)),Tr/(4^)), m = Q,.. 



where the weight functions, the nodal points, and the weights are not the same as for 
the continuous method. 

4.3. The multi-adaptive mcG(q)-mdG(g) method. The discussion above 
for the two methods extends naturally to using different methods for different com- 
ponents. Some of the components could therefore be solved for using the mcG(g) 
method, while for others we use the mdG(g) method. We can even change methods 
between different intervals. 

Although the formulation thus includes adaptive orders and methods, as well as 
adaptive time-steps, our focus will be mainly on adaptive time-steps. 

4.4. Choosing basis functions and quadrature. What remains in order to 
implement the two methods specified by (|4. 6|) and (|4. 1 1[) is to choose basis functions 
and quadrature. For simplicity and efficiency reasons, it is desirable to let the nodal 
points for the nodal basis coincide with the quadrature points. It turns out that for 
both methods, the mcG^) and the mdG(g) methods, this is possible to achieve in 
a natural way. We thus choose the q + 1 Lobatto quadrature points for the mcG(g) 
method, i.e., the zeros of xP q (x) — P q -i(x), where P q is the gth-order Legendre poly- 
nomial on the interval; for the mdG(g) method, we choose the Radau quadrature 
points, i.e., the zeros of P q {x) + P q +i(x) on the interval (with time reversed so as to 
include the right end-point). See [28] for a detailed discussion on this subject. The 
resulting discrete schemes are related to the implicit Runge-Kutta methods referred 
to as Lobatto and Radau methods; see, e.g., [3]. 

5. Basic properties of the mult i- adaptive Galerkin methods. In this sec- 
tion we examine some basic properties of the multi-adaptive methods, including order, 
energy conservation, and monotonicity. 

5.1. Order. The standard cG(q) and dG(q) methods are of order 2q and 2q + 1, 
respectively. The corresponding properties hold for the multi- adaptive methods, i.e., 
mcG(g) is of order 2q and mdG(g) is of order 2q + 1, assuming that the exact solution 
u is smooth. We examine this in more detail in subsequent papers. 

5.2. Energy conservation for mcG(q). The standard cG(q) method is 
energy-conserving for Hamiltonian systems. We now prove that also the mcG(q) 
method has this property, with the natural restriction that we should use the same 
time-steps for every pair of positions and velocities. We consider a Hamiltonian sys- 
tem, 




or, applying nodal quadrature, 



n=0 



(5.1) 



x = -V x P(x) 
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on (0,T] with x(t) G R N , together with initial conditions for x and x. Here x, is the 
acceleration, which by Newton's second law is balanced by the force F(x) = —V x P{x) 
for some potential field P. With u — x and v = x we rewrite (|5.ip as 



(5.2) 



ii 




V 




' fu(v) ' 


V 




F{u) 




_ fv(u) _ 



f(u,v). 



The total energy E(t) is the sum of the kinetic energy K(t) and the potential energy 
P(x(t)), 



(5.3) 
with 

(5.4) 



E(t) = K(t) + P(x(t)), 



^(*) = ^l|i(*)H 2 = ^ll«(*)ll a . 



Multiplying (|5.1I) with x it is easy to see that energy is conserved for the continuous 
problem, i.e., E{t) = E(0) for all t £ [0, T}. We now prove the corresponding property 
for the discrete mcG(g) solution of (|5.2[) . 

Theorem 5.1. The multi- adaptive continuous Galerkin method conserves energy 
in the following sense: Let (U,V) be the mcG(q) solution to (|5.2[) defined by (|4.3p . 
Assume that the same time-steps are used for every pair of positions and corresponding 
velocities. Then at every synchronized time-level i, such as, e.g., T , we have 



(5.5) 



K{t)+P(t)=K{0) + P{0), 



with Kit) = ±\\V{t)\\ 2 and P{t) = P{U{t)). 

Proof. If every pair of positions and velocities have the same time-step sequence, 
then we may choose V as a test function in the equations for (7, to get 



(U, V) dt= / (V, V) dt 



1 d_ 
dt 



\\V\\ 2 dt = K{t)-K(0). 



Similarly, U may be chosen as a test function in the equations for V to get 



(V, U) dt 



-VP(U)U dt = - f j f P{U) dt = -(P(t) - P(0)), 



and thus K(t) + P(t) = K(0) + P(0). □ 

Remark 5.1. Energy conservation requires exact integration of the right-hand 
side / (or at least that J*(U, V) dt + (P(t) - P(0)) = 0) but can also be obtained in 
the case of quadrature; see [23] . 

5.3. Monotonicity. We shall prove that the mdG(q) method is P-stable (see 

0). 

Theorem 5.2. Let U and V be the mdG(g) solutions of (jl.ip with initial data 
U(0~) and V(0~), respectively, defined by (|4.7p on the same partition. If the right- 
hand side f is monotone, i.e., 

(5.6) (f(u,-)-f(v,-),u-v)<0 Vu,v€R N , 
then, at every synchronized time-level i, such as, e.g., T, we have 

(5.7) \\U(i-)-V(i-)\\<\\U(0-)-V(0-)\\. 
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Proof. Choosing the test function as v — W = U — V in (I4.7[) for U and V, 
summing over the local intervals, and subtracting the two equations, we have 



[W, 



id— 1 



/ WiWi dt 


-I 


'hi 





Noting that 

[W i ]i <j -iW£ j _ 1 + J^.WiWi dt 



(f(U,-)-f(V,-),U-V) dt<0. 



Wj-i + 5(TO a -TO-i) a ). 



we get 



so that 



\w(o-)\\ 



-\\W(T- 



< 



Ei 



»>i— i 



WW; < 0, 



l|w(r-)|| < ||w(o-)||. 

The proof is completed noting that the same analysis applies with T replaced by any 
other synchronized time-level i. □ 

Remark 5.2. The proof extends to the fully discrete scheme, using the positivity 
of the quadrature weights. 

6. A posteriori error analysis. In this section we prove a posteriori error 
estimates for the multi-adaptive Galerkin methods, including quadrature and discrete 
solution errors. Following the procedure outlined in the introduction, we first define 
the dual linearized problem and then derive a representation formula for the error in 
terms of the dual and the residual. 

6.1. The dual problem. The dual problem comes in two different forms: a 
continuous and a discrete. For the a posteriori error analysis of this section, we will 
make use of the continuous dual. The discrete dual problem is used to prove a priori 
error estimates. 

To set up the continuous dual problem, we define, for given functions v\(t) and 
v 2 (t), 



(6.1) 



J*(v 1 (t),v 2 (t),t) 



'9/ 
du 



(svi(t) + (l-s)v2(t),t) ds 



where * denotes the transpose, and we note that 



J(vi, v 2 , -)(vi- v 2 ) = Jg §£(svi + (1 - s)v 2 , ■) ds (vi - v 2 ) 



(6 ' 2) = Jo'f (sv 1 + (l-s)v 2 ,-)ds = f(v 1 ,-)-f(v 2 ,-). 

The continuous dual problem is then defined as the following system of ODEs: 

J*(u,U,-)<p + g on[0,T), 



(6.3) 



-<P 



with data ipx and right-hand side g. Choosing the data and right-hand side appropri- 
ately, we obtain error estimates for different quantities of the computed solution. We 
shall assume below that the dual solution has q continuous derivatives (<^ ? ^ € C(Iij) 
locally on interval Zy) for the continuous method and q + 1 continuous derivatives 
((Pi +1 ^ € C(Iij) locally on interval ly) for the discontinuous method. 
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6.2. Error representation. The basis for the error analysis is the following 
error representation, expressing the error of an approximate solution U (t) in terms of 
the residual R(U, t) via the dual solution ip(t). We stress that the result of the theorem 
is valid for any piecewise polynomial approximation of the solution to the initial value 
problem (jl.ll) and thus in particular the mcG(q) and mdG(q) approximations. 

Theorem 6.1. Let U be a piecewise polynomial approximation of the exact so- 
lution u of (|l.lj) . and let ip be the solution to (|6.3|) with right-hand side g(t) and 
initial data ifx> and define the residual of the approximate solution U as R(U,t) = 
U(t) — f{U{t),t), defined on the open intervals of the partitions Ujiy as 

Ri{U,t) = Ui(t) - fi(U{t),t), t g (fc M -i,fc y ), 

j = 1, . . . , Mi, i = 1, . . . , N. Assume also that U is right- continuous at T . Then the 

error e = U — u satisfies 

(6.4) 

r T N Mi 



,(e) = (e(T) 



,tpr)+[ {e,g)dt = YV^ f R, 
J ° i=l 4=1 Y Jl « 



(U,')<fi dt+ [Uijij-tipifcj-i) 



Proof. By the definition of the dual problem, we have using (|6.2p 



J^{e,g) dt = J T (e, -tp - J* (it, U, -)ip) dt 

= Eij ~ e m dt + Jo (~~ J ( u ' u > ') e ' v) dt 

= En S hj dt + f h . (Mu, ■) - mu, dt. 

Integrating by parts, we get 



hi 



-enpi dt = e i (tl j _ i )tp(t itj ^ 1 ) - ei(ty)v(*ij) + / ^ifi dt, 

hi 



so that 



Eij I hj ~ e m dt = EyNij-l^fe-l) - ( e ( T )> <Pt) + J T (e, <p) dt 
= 12ij[ U i]i,3-l ( Pi{ t iJ-l) - ( e ( T )^T) + Jo(e,lf) dt. 



Thus, with Lp Tt g(e) — (e(T),<px) + Sq(&,9) dt, we have 



( e ) = Eij Liei + •) - /i(^, ^ft * + [f/jlij-i^^j-i) 



= T 



which completes the proof. □ 

We now apply this theorem to represent the error in various norms. As before, 
we let || • || denote the Euclidean norm on 1^ and define IMIliqo.tj.r™) = Jq IMI dt. 

Corollary 6.2. Ifip T = e(T)/||e(T)|| and g = 0, then 



(6.5) 



N Mi 

ht)\\= EE 
i=i J= i 



Ri(U, -)cpi dt + [Uilij-itfifcj-!) 
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Corollary 6.3. Ifcp T = Q and g(t) = e(t)/\\e{t)\\, then 



N Mi 

(6-6) IN|ii([ 0) T],H*0 = X]5Z 

t=l 3=1 



Ri(U,-)ipi dt + [Ui]i t j-i(pi(tij-i) 



6.3. Galerkin errors. To obtain expressions for the Galerkin errors, i.e., the 
errors of the mcG(q) or mdG(<7) approximations, assuming exact quadrature and exact 
solution of the discrete equations, we use two ingredients: the error representation 
of Theorem 16.11 and the Galerkin orthogonalities, (|4.4j) and (14.91) . We first prove the 
following interpolation estimate. 

Lemma 6.4. If f e C q+1 ([a, &]), then there is a constant C q , depending only on 
q, such that 



(6.7) 



1 



\f{x)-^f{x)\<C q k« +x - I \f (q+1) (y)\dy Vie [a, 6], 



where TT^f(x) is the qth-order Taylor expansion of f around Xq = (a + b)/2, ft = b — a, 
andC q = l/(2"q\). 

Proof. Using Taylor's formula with the remainder in integral form, we have 



\f{x)-^f{x)\ = 



A f f iq+1) (y)(y-xo) iq) dy 



□ 



Note that since we allow the polynomial degree to change between different com- 
ponents and between different intervals, the interpolation constant will change in the 
same way. We thus have C qi = C qi (t) = C qij for t e lij . 

We can now prove a posteriori error estimates for the mcG(g) and mdG(g) meth- 
ods. The estimates come in a number of different versions. We typically use E 2 or £3 
to adaptively determine the time-steps and E or E\ to evaluate the error. The quan- 
tities E4 and E$ may be used for qualitative estimates of error growth. We emphasize 
that all of the estimates derived in Theorems 16.51 and 16.61 below may be of use in an 
actual implementation, ranging from the very sharp estimate Eq containing only local 
quantities to the more robust estimate E$ containing only global quantities. 

Theorem 6.5. The mcG(q) method satisfies the following estimates: 



(e)| = E < Ex < E 2 < E 3 < E 4 



and 

(6.9) 

where 



(6.10) 



<E 2 <E 5 



E {) 
Ei 

E-2 

E?, 
Ei 
E, 



Eili s l qz] max [o,T] {c qi -ikf n} . 



TTklfii) dt 

n k ipi\ dt 
fey] 



SW-Wmaxj , [0 ,t] {C^kfn} , 
S^ 2 \\C q ^k q R(U,-)\\ L 2 (RKmT]) , 
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with C q as in Lemma WM ki(t) = kij, ri(t) — r^, and sf*\t) — s,^ !J ' for t 6 1^, 

s [qi - 



(6.11) = / T |^| dt, ^ = jf\\<pM\\dt, 

S M,a = (/ r ||^fdt) 1/a , 

and where irkip is any test space approximation of the dual solution ip. Expressions 
such as C q -\k q R are defined componentwise, i.e., (C q -ik q R(U, •)){ — C qij -\kf? Ri(U, ■ 
for t G Iij . 

Proof. Using the error representation of Theorem 16.11 and the Galerkin orthogo- 
nality (|4.4jl . noting that the jump terms disappear since U is continuous, we have 



JV Mi 



J2J2f j^oOi-w) dt 

i=l 7=1 



= E , 



3- 

whcre i^^ip is any test space approximation of (p. By the triangle inequality, we have 

N Mi 

E <Y,J2 \Ri{U,')(<Pi-*m)\dt = Ei- 

i=l j=l J Ii i 

Choosing Kkfi as in Lemma 16.41 on every interval /y, we have 

Ex ■ y^c, f \iu(u,.)\ « -1 f \^ ] \ dt 

Continuing, we have 



E2 < Etinmx^C^fcfrJE^Mr 
= E*=i max [o,T] {C 3i _ifef r»} J T \(p\ qi) \ dt 

= Eti max io,T] {C qi -ik?n} = e 3 , 



and, finally, 



E 3 < max l , [0 , T] {C g ,_ 1 fcfr l }Eti/o T |^ l) l * 
< max M0 , T] {C^Afr*} VF/* dt 
= max^o.T] {C qi -xkfn} VNS^ 1 = E*. 

As an alternative we can use Cauchy's inequality in a different way. Continuing from 
E2, we have 

E2 = ELE^iCW-i^ +1 ^if ] 

= Eti EfA ^-i^/4f 1 /,„ I W 01 ^ 

- EUoC^kriR^-^dt 

= £(C q ^k*\R(U,-)\,sM)dt 

< r T ||c ? _i^i?([/,-)||||sM|| dt 



< 



^\\c q ^R{u,-)f dt) 1 ' 2 (; T || S M|| 2 dt) 1/2 , 
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where \R(U, -)| denotes the vector-valued function with components = = 
\Ri(U, -)|. Noting now that s is the L 2 -projection of \ | onto the piecewise constants 
on the partition, we have 



1/2 



\\s [q] f dt 



< 



1/2 



M 9) H 2 dt 



so that 



\L 



r,g(e)\ < \\C q -ik q R(U, -)\\l 2 (R n xlO^W^ \\l 2 (R n x[0 ,T]) = -^5, 



completing the proof. □ 

The proof of the estimates for the mdG(g) method is obtained similarly. Since 
in the discontinuous method the test functions are on every interval of one degree 
higher order than the test functions in the continuous method, we can choose a better 
interpolant. Thus, in view of Lemma 16.41 we obtain an extra factor kij in the error 
estimates. 

Theorem 6.6. The mdG(g) method satisfies the following estimates: 



(6.12) 

and 

(6.13) 

where 
(6.14) 

E() 
E x 

E-2 

E3 

Ei 

£5 

with 
(6.15) 



\L VT , g (e)\ = E <E 1 <E 2 <E 3 <E i 



\L VT ,g(e)\ <E 2 < E 5 , 



Y,ij Si v Ri{U, '){<p% - n k <Pi) dt + [Ui]i t j-i(tpi(t it j-i) - Wk<Pi(fij-i)) 

Eij Sl tj {Rifti OllVi - TTfeVil dt+ |[Ut]i,3-l||Vi(*t,j-l) - TfcV»(*iJ_l) 



2-,i=i Z^j=i 



r ij S ij 



max [0jT ] <j C qi k\ 



AT max 



»,[0,T] 



Co, 



= S^+^\\C q k^R(U,-)\\ L2(m 



x[0,T])) 



i2i(£/,-) = W,-)l + Trl[^ki-il 



and we otherwise use the notation of Theorem \6.5i 

Proof. As in the proof for the continuous method, we use the error representation 
of Theorem 16.11 and the Galerkin orthogonality (|4.9p to get 



,(e)| 



■_• J li-i 



= E n 



By Lemma l6.4l we obtain 

< Eij C^r +1 (j hj \Ri(U, ■)\dt+ IN^-il) J- ! hj dt 

< Z ij C qij k% + %s%^=E 2 . 
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Continuing now in the same way as for the continuous method, we have E 2 < E 3 < E4 
and E 2 < E 5 . □ 

Remark 6.1. When evaluating the expressions Eq or E\, the interpolant n^ip does 
not have to be chosen as in Lemma 16.41 This is only a convenient way to obtain the 
interpolation constant. In section 16.61 below we discuss a more convenient choice of 
interpolant. 

Remark 6.2. If we replace Jj \Ri\ dt by max/ ;j \Ri\, we may replace C q by a 
smaller constant C' q . The value of the constant thus depends on the specific way the 
residual is measured. 



6.4. Computational errors. The error estimates of Theorems [6.51 and IBTBI are 

based on the Galerkin orthogonalities (|4.4[) and (|4.9[) . If the corresponding discrete 
equations are not solved exactly, there will be an additional contribution to the total 
error. Although Theorem 16. II is still valid, the first steps in Theorems 16.51 and 16.61 are 
not. Focusing on the continuous method, the first step in the proof of Theorem 16.51 
is the subtraction of a test space interpolant. This is possible, since by the Galerkin 
orthogonality we have 

JV Mi . 
i=l 3=1 

for all test space interpolants 7Tfe<£>. If the residual is no longer orthogonal to the test 
space, we add and subtract this term to get to the point where the implications of 
Theorem 16.51 are valid for one of the terms. Assuming now that (p varies slowly on 
each subinterval, we estimate the remaining extra term as follows: 
(6.16) 



E c 



Ei=l E,=l kijl^ij 1 ' 



dt 



Ef=iEnfci3#i3-lK-i 



where ip is a piecewise constant approximation of ip (using, say, the mean values on 
the local intervals), 



3.17) 



m 



Mi 
3 = 1 



\cpi\ dt = S. 



[o] 



is a stability factor, and we define the discrete or computational residual as 



(6.18) 



n c - = - 
1 » k 



Ri{U,-) dt 



h(u,-) dt 



More precise estimates may be used if needed. 

For the mdG(q) method, the situation is similar with the computational residual 
now defined as 



3.19) 



fi(U,-) dt 



Thus, to estimate the computational error, we evaluate the computational resid- 
uals and multiply with the computed stability factors. 
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6.5. Quadrature errors. We now extend our analysis to take into account also 
quadrature errors. We denote integrals evaluated by quadrature with J. Starting 
from the error representation as before, we have for the mcG(q) method 
(6.20) 

L VT , g {e) = J^(R,ip)dt 

= J T (i?, If ~ TTkP) dt + Jq(R, 7T fe <yj) dt 



rp ~T 

= J (R, ip-TTkp) dt + f (R,Tr k (p) dt + 



rp ~T 

J Q (R,TT k if) dt - J Q (R,U k <p) dt 



rp ~T / ~T p\ 

= Jo ( R > V - ^k^) dt + J Q (R, TT k Lp) dt+lj -J j (f(U, •), w k f) dt 

if the quadrature is exact for Uv when v is a test function. The first term of this 
expression was estimated in Theorem 16.51 and the second term is the computational 
error discussed previously (where J denotes that in a real implementation, (|6.18[) is 
evaluated using quadrature). The third term is the quadrature error, which may be 
nonzero even if / is linear, if the time-steps are different for different components. To 
estimate the quadrature error, notice that 



(6.21) 



(Jo - Jo T ) (m -uw) ^ = £y (/ Jy - j hj ) m -y Wi 



dt 



Ey%^|<Eti5l 0I ma Xj -|7e||, 



where {sf ] }f =1 are the same stability factors as in the estimate for the computational 
error and 

(6.22) n % = h\ I /l(C/ ' ') dt ~ f fl{U ' ' ] A 

is the quadrature residual. A similar estimate holds for the mdG(g) method. 

We now make a few comments on how to estimate the quadrature residual. The 
Lobatto quadrature of the mcG(g) method is exact for polynomials of degree less than 
or equal to 2q — 1, and we have an order 2q estimate for J — J in terms of /^ 2<? \ and 
so we make the assumption TZQ oc h 2 ^ 3 . If, instead of using the standard quadrature 
rule over the interval with quadrature residual TZfj , we divide the interval into 2 m 
parts and use the quadrature on every interval, summing up the result, we will get a 
different quadrature residual, namely 

(6.23) K Qm = -C2 m (k/2 m ) 2q+1 = 2 m( -- 2q ^Ck 2q = 2- 2q Tl Qm - 1 , 

k 

where we have dropped the y subindices. Thus, since |7?. Sm | < \7Z Qm — 7Z Qm+1 \ + 
\ftQm+i \ = |7^S m _ TZ Qm + 1 \ + 2- 2q \K Qm l we have the estimate 

(6.24) \K Q ™ | < 1 -±-^\K Q - - TZ Qm+1 \. 

Thus, by computing the integrals at two or more dyadic levels, we may estimate 
quadrature residuals and thus the quadrature error. 

For the mdG(g) method the only difference is that the basic quadrature rule is 
one order better, i.e., instead of 2q we have 2q+l, so that 

(6.25) \n Q -\ < 1 _^_ 1 _ 2q \n Q --n a ^\. 
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6.6. Evaluating Eq. We now present an approach to estimating the quantity 
J (R(U, -),(p — TTk^p) dt by direct evaluation, with ip a computed dual solution and 
■K^Lp a suitably chosen interpolant. In this way we avoid introducing interpolation 
constants and computing derivatives of the dual. Note, however, that although we do 
not explicitly compute any derivatives of the dual, the regularity assumed in section 
16.11 for the dual solution is still implicitly required for the computed quantities to 
make sense. Starting now with 



(6.26) 



En = 



N Mi 

EE 

i=i j=i 



Ri(U, -)((pi - TTkVi) dt 



for the continuous method, we realize that the best possible choice of interpolant, if 
we want to prevent cancellation, is to choose such that Ri{U, -)((pi — i^k^Pi) > 
(or < 0) on every local interval Zy. With such a choice of interpolant, we would have 
(6.27) 

TV Mi 



E G 



N Mi 

EE 

i=i j=i 



Ri(U, -)((pi - n k ipi) dt 



EE- 

»=i i=i 



\Ri(U,-)(ipi-Tr k Vi)\ dt 



with otij — ±1. The following lemmas give us an idea of how to choose the interpolant. 

Lemma 6.7. If, for i = 1, . . . ,N, fi = f t (U(t),t) = fi(Ui(t),t) and fi is linear 
or, alternatively, f — f(U(t),t) is linear and all components have the same time- 
steps and order, then every component Ri(U, •) of the mcG(g) residual is a Legendre 
polynomial of order qij on 1^ , for j = 1, . . . , Mj. 

Proof. On every interval the residual component Ri(U,-) is orthogonal to 
r P qii ~ 1 (Iij). Since the conditions assumed in the statement of the lemma guarantee 
that the residual is a polynomial of degree on every interval Iij, it is clear that 
on every such interval it is the g^th-order Legendre polynomial (or a multiple 
thereof). □ 

Even if the rather strict conditions of this lemma do not hold, we can say some- 
thing similar. The following lemma restates this property in terms of approximations 
of the residual. 

Lemma 6.8. Let R be the local L 2 -projection of the mcG(q) residual R onto 
the trial space, i.e., Ri(U, ■) |/ i . is the L 2 (Iij)-projection onto V qij (hj) of Ri(U, ■) |/ i . , 
j = I,..., Mi, i — 1, . . . , N. Then every Ri(U, ')|i 4;f is a Legendre polynomial of degree 
qij. 

Proof. Since Ri(U, ■) is the L 2 -projection of Ri(U, •) onto V qij (Iij) on Z^, we have 
Ri{U,-)vdt= R l (U,-)vdt = 

for all v € V qi 3~ l {Iij)i so that Ri(U,-) is the g^-th-order Legendre polynomial on 

hj. o 

To prove the corresponding results for the discontinuous method, we first note 
some basic properties of Radau polynomials. 

Lemma 6.9. Let P q be the qth-order Legendre polynomial on [—1,1]. Then the 
qth-order Radau polynomial, Q q (x) = (P q {x) + P q +i(x))/(x + 1), has the following 
property: 



(6.28) 



/ - / Q q {x)(x + l) p dx = 
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for p = 1, . . . , q. Conversely, if f is a polynomial of degree q on [— 1, 1] and has the 
property (|6.28[) . i.e., f(x)(x + l) p dx = for p = 1, . . . ,q, then f is a Radau 
polynomial. 

Proof. We can write the gth-order Legendre polynomial on [—1,1] as P q (x) — 
-^2vD q {{x 2 — l) q ). Thus, integrating by parts, we have 

j = ji^ p,M+P|+iW ( i + 1 )p dx 

= S-i DQ ^ x2 ~ *)* + x ( x% ~ l ) q )& + !) P_1 dx 

= ilk I-i Dq (( x + !)^ 2 - !)«)(* + dx 

= f-i D q - p ((x + l)(x 2 - lY)DP(x + If- 1 dx = 0, 

since D ((x + l)(x 2 — l) q ) is zero at —1 and 1 for / < q. Assume now that / is 
a polynomial of degree q on [—1, 1] with the property (|6.28p . Since {(x + l) p }p =1 
are linearly independent on [—1, 1] and orthogonal to the Radau polynomial Q q , 
{Q q (x), (x + 1), (x + l) 2 , ...,{x+ l) q } form a basis for V>([-1, 1]). If then / is 
orthogonal to the subspace spanned by {(x + l) p }p =1 , we must have / = cQ q for some 
constant c, and the proof is complete. □ 

Lemma 6.10. If, for i = 1, . . . , N, fi = fi{U(t),t) = fi(Ui(t),t) and fa is linear 
or, alternatively, f = f(U(t),t) is linear and all components have the same time- 
steps and order, then every component Ri{U,-) of the mdG(q) residual is a Radau 
polynomial of order q^ on Iij for j = 1, . . . , Mj. 

Proof. Note first that by assumption the residual Ri(U, •) is a polynomial of degree 
qij on Iij . By the Galerkin orthogonality, we have 

= jT R t (U, > dt + [Uilij-Mttj-i) w e 

which holds especially for v{t) = (t — Uj-i) p with p — l,...,q, for which the jump 
terms disappear. Rescaling to [—1,1], it follows from Lemma [6.91 that the residual 
Ri(U,-) must be a Radau polynomial on ly. □ 

Also for the discontinuous method there is a reformulation in terms of approxi- 
mations of the residual. 

Lemma 6.11. Let R be the local L 2 -projection of the mdG(q) residual R onto 
the trial space, i.e., Ri(U, >)|j 4 . is the L 2 (Iij) -projection onto V qij (hj) of Ri(U, 
j = 1, . . . , Mi, i — 1, . . . , N . Then every Ri(U, * s a Ro-dau polynomial of degree 
qij. 

Proof. Since Ri(U, •) is the L 2 -projection of Ri(U, •) onto V qij (Iij) on 1^, it follows 
from the Galerkin orthogonality that 

/ R l (U,-)vdt= Ri(U,-)vdt = 
Jin Jin 

for any v(t) = (t — ti.j-i) p with 1 < p < q. From Lemma [6.91 it then follows that 
Ri(U, ■) is a Radau polynomial on Jy. □ 

We thus know that the mcG(g) residuals are (in the sense of Lemma l6.8[) Legendre 
polynomials on the local intervals and that the mdG(g) residuals are (in the sense of 
Lemma |6. lip Radau polynomials. This is illustrated in Figure RTT1 

From this information about the residual, we now choose the interpolant. Assume 
that the polynomial order of the method on some interval is q for the continuous 
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-2 



-3 - 

5 5.1 5.2 5.3 5.4 , 5.6 5.7 5.8 5.9 6 




-6 



-8 



5 5.1 5.2 5.3 5.4 ^ 5.6 5.7 5.8 5.9 6 

Fig. 6.1. The Legendre-polynomial residual of the mcG(g) method (left) and the Radau- 
polynomial residual of the mdG(q) method (right), for polynomials of degree five, i.e., methods 
of order 10 and 11, respectively. 



method. Then the dual should be interpolated by a polynomial of degree q—l, i.e., we 
have freedom to interpolate at exactly q points. Since a gth-order Legendre polynomial 
has q zeros on the interval, we may choose to interpolate the dual exactly at those 
points where the residual is zero. This means that if the dual can be approximated 
well enough by a polynomial of degree q, the product Ri(U, — itkVi) does not 
change sign on the interval. 

For the discontinuous method, we should interpolate the dual with a polynomial 
of degree q, i.e., we have freedom to interpolate at exactly q + 1 points. To get rid 
of the jump terms that are present in the error representation for the discontinuous 
method, we want to interpolate the dual at the beginning of every interval. This 
leaves q degrees of freedom. We then choose to interpolate the dual at the q points 
within the interval where the Radau polynomial is zero. 

As a result, we may choose the interpolant in such a way that we have 
(6.29) 



\ L <p T ,g( e )\ = 



ij 



Ri(U, -)(ipi - ir k ipi) dt 



ij 



a. 



\Ri(U, - TfeVi)! dt, 
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with ctij = ±1, for both the mcG(q) method and the mdG(g) method (but the 
interpolants are different). Notice that the jump terms for the discontinuous method 
have disappeared. 

There is now a simple way to compute the integrals Jj Ri(U, -)(ifi — T^kVi) dt. 
Since the integrands are, in principle, products of two polynomials for which we know 
the positions of the zeros, the product is a polynomial with known properties. There 
are then constants C q (which can be computed numerically), depending on the order 
and the method, such that 



(6.30) / \Ri(U, - W )| dt = C q MRi{U,t^)\\<pi(tij) ~ 

Finally, note that there are "computational" counterparts also for the estimates 
of type E$ in Theorems 16.51 and 16.61 namely 

\L VTtg {e)\ < EijX/y \Ri{U,-)\\<Pi-irkVi\ dt 

(6.31) = Y lii C' m k%\R i {U i %)\J Iij - i ^\<p i -* kt p i \dA 

< Eli ^ma*, : w C'^kflRiiU,^)], 

with Si = Jq jlri\ipi — Ttk<Pi\ dt for the continuous method and similarly for the 
discontinuous method. 

6.7. The total error. The total error is composed of three parts — the Galerkin 
error, Eq, the computational error, Ec and the quadrature error, Eq: 

(6.32) \L VT , g {e)\ <E G + E C + E Q . 

As an example, choosing estimate E3 of Theorems 16.51 and 16.61 we have the following 
(approximate) error estimate for the mcG(g) method: 



(6.33) \L VT , g (e)\<J2 



N 

Sl gi] max{C„ i - 1 kf i n} + s\ 0] max Iftfl + s\ 0] maxl^l 

[0,T] [0,T] [0,T] 



for the mdG(g) method we have 
(6.34) 

\L VT , g {e)\ < y S\ qi+1] m ax \C qi kf +V;} + ^ [0] max \Kf | + S l ° ] max \K 

* 9 •f-' fO.Tl I J 1 \0,T] 1 [0,T1 



i=l 



These estimates containing Galerkin errors, computational errors, and quadrature 
errors also include numerical round-off errors (included in the computational error). 
Modelling errors could also be similarly accounted for since these are closely related 
to quadrature errors, in that both errors can be seen as arising from integrating the 
wrong right-hand side. 

The true global error may thus be estimated in terms of computable stability 
factors and residuals. We expect the estimate for the Galerkin error, Eq, to be quite 
sharp, while Ec and Eq may be less sharp. Even sharper estimates are obtained 
using estimates Eq, E%, or E2 of Theorems 16.51 and [ 
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6.8. An a posteriori error estimate for the dual. We conclude this section 
by proving a computable a posteriori error estimate for the dual problem. To compute 
the stability factors used in the error estimates presented above, we solve the dual 
problem numerically, and we thus face the problem of estimating the error in the 
stability factors. 

To demonstrate how relative errors of stability factors can be estimated using the 
same technique as above, we compute the relative error for the stability factor Sp(T), 
defined as 



3.35) 



S, 



V (T) = sup 

\MT)\\ = 1 



m 



dt 



for a computed approximation $ of the dual solution tp. 

To estimate the relative error of the stability factor, we use the error representa- 
tion of Theorem 16. II to represent the i 1 ([0, T], E JV )-error of $ in terms of the residual 
of $ and the dual of the dual, w. In [25J we prove the following lemma, from which 
the estimate follows. 

Lemma 6.12. Let (p be the dual solution with stability factor S v (t), i.e., with 
data \\ip(t)\\ = 1 specified at time t, and let uj be the dual of the dual. We then have 
the following estimate: 



(6.36) 



||w(t)|| < S v (T-t) Vie[0,T]. 



Theorem 6.13. Let $ fee a continuous approximation of the dual solution with 
residual and assume that S v (t)/ S V (T) is bounded by C on [0,T]. Then the 
following estimate holds for the relative error of the stability factor S<& (T) : 



\S*(T) - S V (T)\/S V (T) < C / dt, 



3.37) 



and for many problems we may take C = 1. 

Proof. By Corollarv l6.3[ we have an expression for the i 1 ([0,T],]R Ar )-error of the 
dual, so that 



f T (||<i>||-|MI)^ <fi\\$-v\\dt 



\S*(T) - S V (T)\ = 

(6.38) 

< ; J ||jz*|||Kr-.)|| dt 

With C defined as above it now follows by Lemma 16.121 that 

l-T 



* - </>IUi([0,T],K») =Jo ( R <S>^( T ~ ■)) dt 
T i 



|S«(T) - S V (T)\ < C [ \\R*\\ dt S V (T), 
Jo 



and the proof is complete. □ 

Remark 6.3. We also have to take into account quadrature errors when evaluating 
(|6.35[) . This can be done in many ways; see, e.g., [3]. 

Appendix A. Derivation of the methods. 

This section contains some details left out of the discussion of section 0J 
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A.l. The mcG(q) method. To rewrite the local problem in a more explicit 
form, let {s n }^ =0 be a set of nodal points on [0,1], with s = and s q = 1. A good 
choice for the cG(q) method is the Lobatto points of [0, 1]. Now, let Ty be the linear 
mapping from the interval to (0, 1], defined by 

(A.l) Tij (t)- 



and let {A$}^_ be the {sn}' =0 Lagrange basis functions for V q ([0, 1]) on [0, 1], i.e. 

(A 2) A M (s) = ( s ~ s o) ' ' ' ( s " fn-lKj - gn+l) ■ ■ • (f - gg) 

(fin - so) •• • (s n — s n _i)(s n - s„+i) •• • (s n - s 9 ) 



We can then express E/j on Zjj in the form 

(a.3) m) = J2^^ ] (n j (t)), 



and choosing the A m 1 ' as test functions we can formulate the local problem (|4.3p as 
follows: Find {fy n }n=o' witn &i0 = &,j-i,9< 3 -i i such that for m = 0, . . . , %j - 1 

(A - 4) ? 

/ N? y W*))l / fi(U(t),t)>&- 1] (ni{t)) dt. 

J I*i n=0 dt 1 J Jl » 

To simplify the notation, we drop the ^ subindices and assume that the time-interval 
is [0,k], keeping in mind that, although not visible, all other components are present 
in /. We thus seek to determine the coefficients {£ n }n=i with £ given, such that for 
m = 1 , . . . , q we have 

(A.5) £>i ^ ^ ] (r(t))\ [ ^(T(t)) dt= f f\t--Kr{t)) dt, 

or simply 

(A.6) = &m, 

n=l 

where 

(A.7) a® n = I' X^(t)X [ r-Ht)dt 



and 

(A.8) b m = I /A^Ii ] (t(*)) dt - a mQ £ . 



We explicitly compute the inverse ^4.^ — (o>mn ) of the matrix A® = (a m ] „). Thus, 
switching back to the full notation, we get 



(A.9) £ym = -£„ 5 ™« a »° + / w m j] (nj (*)) /i(I7(t),t) dt, 
n=l 



m = 1,.. 
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where the weight functions {wm}m—i are given by 
(A.10) w^=J2^n^ ] , m = l,...,q. 

n=l 

Following Lemma I A. 1 1 below, this relation may be somewhat simplified. 
Lemma A.l. For the mcG(g) method, we have 

i 

E S ™« a "° = _L 

n=l 

Proof. Assume the interval to be [0, 1]. The value is independent of / so we may 
take / = 0. We thus want to prove that if / = 0, then £„ = £o for n = 1, . . . , q, 
i.e., U = Uq on [0, 1] since {a!?'}* =0 is a nodal basis for P q {[0, 1]). Going back to the 
Galerkin orthogonality (|4.4p . this amounts to showing that if 

Uvdt = Vue7> 9_1 ([0,l]), 



with U £ V q ([0, 1]), then U is constant on [0,1]. This follows by taking v = U. 
Thus, = £o for n = 1, . . . , q, so that the value of J2n=i "mUnO must be —1. This 
completes the proof. □ 

The mcG(g) method thus reads as follows: For every local interval Iij, find 
{&jn}n=0> with 6y = &,3-i,g 4 ,j_i) such that 

(A.ll) £ ijm = Ci i0 + f ( nj {t)) fi(U(t),t) dt, m = 1, ... , qij , 

J la 

for certain weight functions {uin^} q n=1 C 7' 9_1 (0, 1), and where the initial condition 
is specified by &00 = «i(0) for i = 1, . . . , N. 

The weight functions may be computed analytically for small q, and for general 
q they are easy to compute numerically. 

A. 2. The mdG(q) method. We now make the same ansatz as for the contin- 
uous method, 

(A.12) U i (t) = ^2^ jn X^\r ij (t)), 

n=0 

where the difference is that we now have q + 1 degrees of freedom on every interval, 
since we no longer have the continuity requirement for the trial functions. We make 
the assumption that the nodal points for the nodal basis functions are chosen so that 

(A.13) s q = 1, 

i.e., the end-point of every subinterval is a nodal point for the basis functions. 

With this ansatz, we get the following set of equations for determining {£ijn}nLo'- 



(A.14) 

E X n A (°) - £«o J x t ] (0) + / E Tt x t s] (ns (*)) ^Hn 3 (t)) dt 

n=0 J Jl H n=0 



(A.15) = / MU(t),t)X l ^ ] (r tJ (t)) dt 
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for m = 0, . . . , qij, where we use to denote > i- e -j the value at the right 

end-point of the previous interval. To simplify the notation, we drop the subindices 
ij again and rewrite to [0, k]. We thus seek to determine the coefficients {£, n }n=o sucn 
that for m = 0, . . . , q we have 

(A.16) 

(l>AM(0)-£- ] A$(0) + £>1 / fc AM(r(t))AM(r(t)) dt = f f^{r(t)) dt, 
or simply 

(A.17) = &m, 



n=0 



where 



(A.18) a{l = /' Ak 9] (i)Afe 1 (t) d* + AM (0)A$ (0) 

and 

(A.19) 6M = / /AM(r(t)) dt + &*®(0). 

Jo 

Now, let AM be the (q+l) x (qr+1) matrix A^ = (a$„) with inverse AM = (a$„). 
Then, switching back to the full notation, we have 
(A.20) 

9 r 

Ziim = £yo E «mUi? ] (0) + / /*(!/(*),*) dt, m = 0, . . 

n=0 ^ 7 »J 

where the weight functions {u>M}* =0 are given by 

(A.21) ^E^W? 1 ' "» = 0,...,«. 



n=0 



As for the continuous method, this may be somewhat simplified. 
Lemma A. 2. For the mdG(q) method, we have 



£aM„AM(0) = l. 



n=0 

Proof. As in the proof for the mcG(q) method, assume that the interval is [0, 1]. 
Since the value of the expression is independent of / we can take / = 0. We thus want 
to prove that if / = 0, then the solution U is constant. By the Galerkin orthogonality, 
we have 

[U] v(0)+ [ Uvdt = Q VveV q (0,l), 
Jo 

with U e V q (0, 1). Taking v = U- U(0~), we have 

o = mo) 2 + foU(u-u(o-))dt = ([u} ) 2 + y o 1 £(u-u(o-)fdt 

= i(C/(0+)-[/(0-)) 2 + i([/(l)-C/(0-)) 2 , 
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so that [U] = 0. Now take v = U. This gives Q{U) 2 dt = 0. Since then both 
[U]o — and U = on [0, 1], U is constant and equal to U(0~), and the proof is 
complete. □ 

The mdG(g) method thus reads as follows: For every local interval ly, find 
{£ijn}n=o> sucn that for m — 0, . . . , g,j we have 

(A.22) (U im = + I (nj (*)) fi(U(t), t) dt 

for certain weight functions {wn'}* =0 C V q (0, 1). 
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